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Abstract 

We introduce new "convolution spline" methods for the temporal approximation of time domain boundary 
integral equations (TDBIEs). They share some properties of convolution quadrature, but instead of being 
based on an underlying ODE solver are explicitly constructed in terms of basis functions which have 
compact support, and hence give sparse system matrices. They are derived and analysed for convolution 
Volterra integral equations (VIEs): at time step t„ = nh the VIE solution is approximated in a "backward 
time" manner in terms of basis functions (j>j by u{tn — f) ~ X]j=o ''^ri-j4ij{t/h) for t G [0, t„]. We show that 
using isogeometric B-splines of degree m > 1 on [0, co) in this framework gives a second order accurate 
scheme, but cubic splines with the "parabolic runout" conditions at t = are fourth order accurate. We 
establish a methodology for the stability analysis of VIEs and demonstrate that the new methods are stable 
for non-smooth kernels which are related to convergence analysis for TDBIEs, including the case of a Bessel 
function kernel oscillating at frequency 0{l/h). Numerical results for VIEs and for TDBIE problems on 
both open and closed surfaces confirm the theoretical predictions. 

Keywords: Convolution quadrature, Volterra integral equations, time dependent boundary integral equa- 
tions 

AMS(MOS) subject classification: 65R20, 65M12 



1 Introduction 

We derive and test new "convolution spline" time-stepping methods for the time-dependent boundary integral 
equation (TDBIE) 

1 f u(x' ,t—\x' — x\) , , , , , , , 

' ^ ' ' dx' = a{x,t) fora;er, i>0 (1.1) 



4tt Jy \x' — x\ 

for u. It is the single layer potential equation for acoustic scattering from the surface F C M"^ with zero 
Dirichlet boundary conditions and (known) incident field —a{x,t), and is equivalent to 



[ [ k{x'-x,t-t')u{x',t') dx' dt' = a 
Jo Jr 



{x,t) for k{z,t) 



Sit-\z\) 



Two key properties of the new methods are: (i) they are easy to implement, and (ii) they are computationally 
efficient (because the resulting system matrices are sparse). They share some properties with convolution 
quadrature (CQ), and we follow Lubich (1988) in deriving them as approximation schemes for the convolution- 
kernel Volterra integral equation (VIE) 

t 

K{t')u{t-t')dt' ^a{t), te[0,T]. (1.2) 







Convergence and stability properties of the new schemes are analysed for VIEs, but the focus of the paper is 
not on deriving new methods for VIEs (of which there are already very many), but on using the insight gained 
from VIEs to derive new methods which have good properties for TDBIEs. 
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1.1 Properties of TDBIE approximations 



Designing a good approximation scheme for the TDBIE is nontrivial; challenges include ensuring that 

it is numerically stable, it is not prohibitively hard to implement for a given scattering surface F, and its 
computational complexity is not infeasibly high. There are good survey articles on approximation methods 
for TDBIEs by Ha-Duong (2003) and Costabel (2004), and wc briefly summarise the pros and cons of some of 
the main approaches. 

Bamberger and Ha Duong (1986) proved that a full Galerkin approximation of in time and space is 

stable and convergent for smooth, closed F (this was extended to the case of open, flat F in Ha-Duong (1990)), 
but the stability of the method relics on all the integrals being evaluated very accurately (the key insight 
on how to do this was provided by Terrasse (1993)). In practice this involves converting five dimensional 
volume integrals over irregular (non-polygonal) sub-regions of F x F x [0, T] to surface integrals which arc then 
evaluated using high precision quadrature, and is extremely complicated to successfully implement in practice, 
even for relatively simple F. Collocation schemes for (|1.1|) are far more straightforward to implement, but there 
is little rigorous convergence analysis for them, and numerical instability is often an issue. Methods which 
use a Galerkin approximation in space and CQ in time have obvious attractions: they are based on rigorous 
theoretical analysis (Lubich 1994) and are relatively straightforward to implement. They are also inherently 
far more stable than those which use Galerkin or collocation time approximations (Lubich (1994) showed that 
the CQ method remains stable when the inner product integrals are approximated), but unfortunately the 
disadvantage this time is higher computational complexity. 

All three approaches approximate (jl.ip as a convolution sum of the form X^m^o^™— " ~ s"' which is 
rearranged to give the time-stepping scheme 

n 

qO ^^n_J2 Q" [/"-™ (1.3) 
m— 1 

for U_" e M^'S , the representation of the spatial approximation of u at or near time t" — nh, where the 
right-hand side vector a" is derived from a{x,t). In the case of both Galerkin and collocation approxima- 
tions the matrices Q™ € r^sxJVs sparse - the number of nonzero elements per row of matrix Q'" is 
0{mhi{m, Ng^^}). In particular this means that (|1.3p can be solved in 0{Ng^^) operations once the right- 
hand side is known, and the overall computational complexity to obtain the approximate solution up to time 
Nt h is 0(min{iV|, A^5, N^^^}) operations. For these hyperbolic problems it is usual to use a timestep h 
commensurate with the side Aa; of a typical space mesh element, and in this case Nt ~ N and Ns ~ N'^ for 
N ^ 1/A.T, and the total computational complexity is 0{N^). Although this compares somewhat unfavourably 
with the 0(N'^) computational complexity of a finite difference or finite clement approximation of the PDE 
formulation of the acoustic wave equation in R^, the plane wave "fast" methods developed by Michielssen 
and co-workers (Ergin, Shanker & Michielssen 1999, Ergin, Shankcr & Michielssen 2000, Lu, Wang, Ergin & 
Michielssen 2000) reduces the complexity to 0{N'^ log^iV). 

Using CQ in time results in a solution algorithm (|1.3p in which the matrices Q™ are dense, because the 
underlying basis functions are global (see e.g. Sec. [2] below, Hackbusch, Kress & Sauter (2009), or Banjai 
(2010) for more details), which increases the computational complexity to 0{Ng N"^). The issue is not solving 
(jl.3p for f/" (which can typically be done efficiently by approximating Q'^ appropriately), but in performing 
the matrix-vector products needed to calculate the right-hand side. Lubich explains that the technique of 
Hairer, Lubich & Schlichte (1985) can be used to reduce the overall complexity to 0{Ng Nt log'^ Nt), i.e. 
0{N^ log^ N). A systematic cut-off strategy to replace small matrix entries by zero is described and carefully 
analysed in Hackbusch et al. (2009), and this reduces the storage costs of the method. This is combined 
with panel clustering by Kress & Sauter (2008) to further reduce the storage costs. CQ methods which 
are based on underlying Runge-Kutta ODE solvers have been developed and analysed for TDBIEs (Banjai 
& Lubich 2011, Banjai, Lubich & Melenk 2011). There are several advantage of these methods over linear 
multistep CQ methods: the basis functions are more highly concentrated (Banjai 2010, Figs 1-2), which makes 
sparsifying the Q™ matrices more straightforward; and higher order accurate methods in time are possible. 
Banjai (2010) uses this approach to develop a practical, parallelizable solution algorithm for (|l.ip which he 
illustrates with a number of realistic large-scale numerical examples. 

1.2 New convolution spline methods 

We describe how to construct a new "convolution spline" method for timestepping (|l.ip which shares some 
properties of CQ, but instead of being based on an underlying ODE solver is explicitly constructed in terms 



of basis functions which have compact support. This means that the Q'" system matrices in (|1.3p for our new 
method have the same sparsity pattern as for the Galerkin or collocation approximations described above, 
giving a TDBIE solution scheme whose overall complexity is 0(rmn{N^ Ns, N^^^}) = 0{N^) operations 
(and which could also be potentially speeded up using "fast" methods). It is also far easier to implement than 
the full space-time Galerkin approach. 

We derive the new approximation as a solution method for the VIE (|1.2p . with u approximated in terms 
of B-spline basis functions in a "backwards time" framework. Our initial approach is to use isogeometric 
B-splines of degree m on [0,00). There can be advantages in using higher order values of m even though the 
formal convergence rate of this scheme for a smooth VIE problem is limited to second order (because it is 
based on quasi- interpolation by the Schoenberg B-spline operator). For example, as noted in Sauter & Veit 
(2011a), using smooth temporal basis functions greatly simplifies approximating the integrals in (jl.ip . We also 
consider cubic B-splines with the "parabolic runout" condition at t = and show that these are fourth order 
accurate. We carefully test out the new methods on ()1.2p . establishing formal convergence, and examining the 
behaviour for kernels which mimic some of the important properties of TDBIE problems, such as discontinuous 
step-function kernels (see e.g. Sauter & Veit (20116)). Another important test problem is obtained from taking 
the spatial Fourier transform of (|l.ip at frequency a; G when F = M^. This is 

[ Jo(ut')u{uj,t-t')dt' = 2a{uj,t), (1.4) 
Jo 

where w = |a;| and Jo is the first kind Bessel function of order zero. As noted in Davies (1994), instabilities 
of approximation schemes for (jl.ip are typically exhibited at the highest spatial frequency which can be 
represented on the mesh. Hence it is important to ensure that any prototype numerical scheme for time- 
stepping (|l.ip is stable for (|1.4p at values of cj = 0{l/h) (assuming h « Ax). 

1.3 Outline 

Section [5] contains an alternative derivation of Lubich's (1988) CQ method for (|1.2p in terms of basis functions 
which have the sum to unity property (j2.12p . The new "convolution spline" approximation of (jl.2p is described 
in SectionOin terms of basis functions which have compact support and are (essentially) all translates, and we 
give sufficient conditions for this approximation to be stable. We consider the case in which the basis functions 
are mth degree isogeometric B-splines on [0,oo) in Section 31 showing how Laplace transform techniques can 
be used to prove the stability of this approximation of (|1.2p for several different test kernels, and demonstrating 
second order convergence for (|1.2p when K and a satisfy 

a e C"^+^[0,r] , K e C'^+^[0,T] , a(0) = and A'(0) = 1 (1.5) 

for suitable d> Q. Under these assumptions, equation (|1.2p possesses a unique solution u E C^lO, T] - e.g. see, 
Brunncr (2004, Theorem 2.1.9). 

In Section [5] we consider a cubic convolution spline basis which is modified near t = to satisfy the "parabolic 
runout" conditions, and show that this gives a far more stable approximation of (jl.2p which is fourth order 
convergent. Numerical tests indicate that this new scheme performs far better than CQ based on BDF2, 
exhibiting fourth order accuracy even for a discontinuous kernel. 

We present numerical test results for TDBIEs in Section [6] which use a Galerkin approximation in space (based 
on triangular piecewise constant elements), and the cubic convolution spline basis in time, for both open and 
closed surfaces F. The TDBIE test problems are similar to those considered in Davies & Duncan (2013) 
which use the convolution-in-time framework with non-polynomial (global) basis functions, but the modified 
B-spline basis functions give a more accurate temporal approximation. 

2 CQ based on linear multistep methods for (11.21) 

We begin by outlining the derivation of the CQ method for (|1.2p from Lubich (1988), in order to show how it 
can be reinterpreted in terms of CQ basis functions. For simplicity we restrict attention to the case for which 
the extension of the solution u by zero to the negative real axis is in C'^{—oo,T] (otherwise the CQ method 
needs to be 'corrected' as described in Lubich (1988, Sec. 3) in order to attain optimal convergence). This is 
guaranteed by requiring 

a(P)(0) = forp = 0:d+l (2.1) 



because u(p)(0) = a'-P+^^O) -J^^Z^ K^p-'^\0) u(^)(0)- We also assume that the Laplace transform J!:(s) of the 
kernel K is sufficiently well-behaved for all the formal manipulations in the next subsection to be rigorous. 
For details see for example (Lubich 1994, Sec. 1) or (Banjai 2010, App). 



2.1 Lubich's CQ method 

We follow Lubich (1988) and substitute the Laplace inversion formula for K{s) into (|1.2p to obtain 

ait) = ^ fK{s)y{t,s)ds, (2.2) 

where 7 is an infinite contour within the region of analyticity of K{s) and y{t, s) = e^* u{t — t') dt' . Treating 
the Laplace variable s as a parameter, y{t) solves the ODE: 

yit) = syit) + u{t), y(0) = 0, (2.3) 

and this is approximated by the fc— step (fc < d) linear multistep method with timestep h 

k k 
j=0 j=0 

where tn = nh, ?/„ sa y{tn) and /„ = sy„ + u{tn). The starting values are y^k = . . .y-i = because of the 
assumption (|2.ip . Multiplying (|2.4p by ^" and summing over n (for ^ G C for which the sum converges) gives 

/^/^\\oo 00 k / ^ 

nr ~ ' ) ^ ^" ^ ^ "(^"^ ^" ' ^^"'^ ^^^^ ^ 5Z / E ft ?'~' 

^ ^ n=0 n=0 j=0 / j=0 

is the symbol of (|2.4p . Hence t/„ is the coefficient of ^" in the expansion of (^^^ — X^fc^o "(^fc) ^'^ ■ 
Substituting ?;„ for y(t„) in (|2.2p shows that a(f„) is approximated by the coefficient of ^" in 

2^z /(^"') '^(■^)^-^E"(^fe)e' = i?('5(e)//i)EH(t.)^'= 



/c=0 /c=0 



using Cauchy's integral formula. Hence, defining the CQ weights qk = qk{h) to be the coefficients in the 
expansion 

CJO 

K{5{0lh) = Y.'i^e (2.5) 

fe=0 

gives the CQ approximation of (|1.2p 

n 

This can be rearranged to give the time-stepping approximate solution m„ sa w(in) 



90 



- (a(t„)-VgjM„_j I forri>l, (2.7) 



since by assumption uq = m(0) = 0. 



2.2 Derivation of CQ in terms of basis functions 

The CQ approximation scheme ()2.6p for the VIE (|1.2p is defined solely in terms of the weights qk- But if CQ is 
used to time-step a TDBIE, then the approximation involves "CQ basis functions" - see e.g. (Hackbusch et al. 
2009, Banjai 2010, Monegato, Scuderi & Stanic 2011). However, we are not aware of a general interpretation 
of CQ approximation schemes for (|1.2p in terms of basis functions. As well as yielding some interesting 
observations, this also gives the framework which we use for the derivation of our "convolution spline" methods 
in Sections mo 



At t = tn := nh (|1.2p can be written as 



a{tn) = / K{t') u{tn - t') dt' , (2.5 



because u{t) = for t < 0. We show below that the standard CQ method is equivalent to approximating u in 
(1^ by 

n 

u(i„ - t') ~ X! ^"-J (^'/^) for > (2.9) 

3=0 

where (f>j are basis functions. This is fundamentally different from standard finite-element type approx- 
imations which have the form u{t) w X^feLo '"fe ''Afe(V^) ~ i-^- the basis function i/jk is always associated with 
the same "unknown" Uk- (Although note that when all the ipk are translates, as happens for some low order 
schemes, then the two approaches can be similar; see Section [3] for more details.) 

Substituting (j2.9|) into (|2.8p and comparing the resulting expression with (|2.6p gives the relationship between 
the standard CQ weights and basis functions: 



qj = / K{t)(t>j{t/h)dt. (2.10) 
Jo 

Comparing this with the standard CQ definition of qj in (|2.5p gives (see Banjai (2010, Eq. (3.1))) 

oo 

e-^(«)*=^0,.(t)e^ (2.11) 
An immediate consequence is that the basis functions satisfy the sum to unity property 

oo 

^0,(t) = l, (2.12) 

provided the underlying multistep ODE solver is consistent, because in this case S{1) — 0. This new observation 
is a crucial property which we use in Section [3] 

2.3 CQ basis functions for LMMs 

Explicit formulae for the 0j(t) based on BDFl-2 have been used for TDBIE approximations (Hackbusch 
et al. 2009, Monegato et al. 2011). The formula for m = 1 is given in Moncgato et al. (2011), and in this case 
<i>j{fy = e~*P /jl, i.e. they arc Erlang functions, used in statistics as probability density functions and satisfy 
<f)j{t) > and 0j(i) dt = 1. The derivation for BDF2 is more complicated, and the explicit formula 

is given in Hackbusch et al. (2009), where Hj is the jth Hermite polynomial. Note that the properties of Hj 
imply that (j}j{t) involves a jth degree polynomial and an exponential in t with no fractional powers of t. 
In principle (j2.1ip can be used directly to find the basis functions (f>j(t) corresponding to any underlying 
linear multistep ODE method for (|2.3p . although this may not be easy in practice. For the trapezoidal rule 
6(0 = 2 (1 - 0/(1 + (Banjai 2010) and is E,°lo (0 = e"'* /(^), where /(C) = exp(4 ^(1 + 0)- 

This gives fiO = E^o f] where 



1 d^ 

J\d^ ^^^^ 



j!(4i).J' dzi 



2=1/(44) 



using the change of variables z = (1 + 0/(4i). It follows from Olver, Lozier, Boisvert & Clark (2010, 
eq. 18.5.6) that fj = (— 1)-' e"''* LJ^(4<) , where L'^{x) is a Laguerre polynomial. The identity Lj^{x) = 
Lj{x) — Lj-i{x) (Gradshteyn & Ryzhik 1994, eq. 8.971-5) gives the trapezoidal rule basis functions 0j(i) = 
(—I)-' {ij{4:t) — £j_i(4t)}, where £j{x) = e^^l"^ ^j{x) is the jth Laguerre function. They are oscillatory, but 
do satisfy (/>j(i) dt = \. These low order basis functions are illustrated in Figure [1] below (see also Banjai 
(2010, Fig. 1) and Monegato et al. (2011, Fig. 4)). 



Scheme 


Initial 

= 0, n > 1 


Recurrence for basis functions 

i> 1 


BDFl 




j</)j(t)-t0,_i(t) =0 


BDF2 


Mt) = e-'"'' 


jcl)jit)^2t<j)jMt) + t^jMt) = o 


BDF3 




j(j)j{t) - 3%_i(t) + 3tq)jMt) - tc^jMt) = 


BDF4 


Mt) = e-''"''' 


iMt) - ^t^,-i{t) + Qt(j)jMt) - 4t0,_3(t) + t(j)jMt) = 


Trap, rule 


Mt) = e-'' 


j0,(t) - Atcb.Mt) + 2(j - l)Miit) + U - 2)0j-2 W = 



Table 1: Recurrence relations for the CQ basis functions. 
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Figure 1: Typical CQ and spline basis functions. See Sections 12.31 and l4.ll for details. 



The direct approach appears intractible for more complicated schemes (even for BDF3), and recurrence rela- 
tions for the basis functions are given in Monegato et al. (2011, Sec. 3.2). They can be compactly derived by 
formally differentiating the generating function (|2.1ip with respect to ^ to get 

oo oo 

i=i 3=0 

and then collecting terms in ^. The initial conditions are (j)n{t) = for n < 0, and the first term of the Taylor 
expansion of (|2.1ip gives (poit) = e^'^'^-'* = e^*"*. Recurrence relations for BDFl-4 and the trapezoidal rule 
are given in Table [1] 



3 Convolution spline approach 

As discussed in Sec. [U basis functions with global support (such as those described above) give rise to dense 
matrices Q™ in the TDBIE scheme (|1.3p . and this has storage and computational cost implications. Here 
we explore the use of compactly supported basis functions, which although not derived via standard CQ, 
nevertheless do fit into the CQ form (j2.9p . We set up a general framework for basis functions for the VIE 
(|1.2p which arc (mainly) translates, and consider specific examples based on B-splines in Sees. I3HS1 This new 
approach gives sparse system matrices when used to time-step TDBIEs, and results are presented in Sec. [SI 

3.1 Construction of a convolution spline scheme for ( 11.21) 

We consider approximations of the form (|2.9p , but where all the basis functions have compact support 
of width 0{h) and almost all are translates of a standard, compactly supported basis function (pm, i-C. 

(f)j{t/h) = (f)„i{t/h + m — j) for j > TO. (3.1) 



The translate property p.ip means that the approximation U{tn — t) fa u{tn — t) has the form 

'""^ /t\ " ft \ 



1 I ' / , "• — I T III \ J 

V ^/ V 

for i > 0; where Vj approximates u{t) for i near (but not necessarily at) tj, and a sum is defined to be zero 
if its upper index is less than its lower index. Substituting the approximation p.2p into the integral equation 
(|1.2p and collocating at each time level as described in Section gives 

n 

J2 Qj Vn-j ^ a{tn) (3.3) 
for n ~ : N where the weights qj are defined by (j2.10p . The unknown coefficients {vj}^^Q are then found by 



time marching as in p.7p . 
3.2 Stability of (IgTSj) 

For TDBIE applications and analysis (see e.g. Davies & Duncan (2004)), we require the scheme p.3p to be 
stable in the following sense, independent of the input function a{t). 

Definition 3.1 (Stability). The scheme i3.3\) is said to be stable when the impulse response sequence {pn} 
defined by 



Po = 1 



1 

~ X! * P''-i forn>l (3.4) 



90 ,=1 



satisfies \pn\ < C for all n such that nh < T , where the constant C is independent of h. 

This is weaker than BIBO (bounded input bounded output) stability in the signal processing literature (see 
e.g. Proakis & Manolakis (2007)), which requires boundedness of the absolute sum X^J^o b"! ^ 
Stability properties of the scheme p.3p can be established by using the Z-transform, defined as follows. 



Definition 3.2. The Z-transform of a sequence {/nj^g function F given by 

F{0^Z{fr.}{0 = J2fnC (3.5) 
n=0 

where ^ € C with \^\ < 1 is such that the sum converges. 
The Z-transform of the scheme (13.31) is 



QiOViO^AiO, (3.6) 

where 

QiO^y,^' K{t)cb,{t/h)dt (3.7) 

and we take a„ — a(t„). Similarly, the Z-transform of (|3.4p is Q{(,)P{0 ~ qo, and so P(^) = qo/Q{0 ■ We 
now state a sufficient condition for stability when Q{^) is a rational polynomial. 

Theorem 3.1 (Root condition for stability). If the Z-transform Q{£,) o/ {g„} is a rational polynomial in ^, 
then the approximation i3.S\} is stable in the sense of Definition \3.I\ if the roots of Q{S^) satisfy the following 
for any constant c > (independent of h): > 1/(1-1- ch) and any with 1/(1 + ch) < |^fc| < 1 are simple. 

Although this result is a variant of the "root condition" familiar (after the change of variable z — 1/^) from 
zero stability analysis of numerical methods for ODEs, it does not appear to have previously been derived or 
used to determine the stability of VIE schemes. Note that simple roots with \£^k\ = 1/(1 + ch) make a bounded 
contribution to Pn as n increases by the standard result 



cT 



iar" = (i + c/i)"<e 

for tn < T, but roots of this size with multiplicity fJ, > 2 grow like n^^^ and hence violate the stability 
definition. 



Verifying the stability condition directly or via the root condition above for a general approximation scheme 
for (|1.2p may be very complicated. But as we show below, schemes with the translate property p.ip can be 
tackled within the framework of Laplace transforms originally introduced for CQ, and this approach gives a 
way to extend the scope of stability analysis to a far broader range of kernel functions. 

Substituting the Laplace inversion formula for K into (|2.10p gives qj = / K{s) $_, (— s/i) ds , where <&, (s) 

is the Laplace transform of (f>j . Hence the approximation scheme p.3p can be written as 

a{tn) = ^ Kis) yn{sh) ds (3.8) 

n 

where yn[sh) = h'^^Vn-j'^j{—sh) . Wc note that ?/„ plays the same role here that the approximate solution 

3=0 

of the ODE (HSl) does in standard CQ. 

The translate property p.ip and the compact support of implies 

^j{~sh) = e'^'**^-™) $m(-s/i) for J > m 

and so 

m 

yn{sh) - e'^yn-ii^sh) = hvn^o{-sh) + h^^Vn-j {^j{~sh) - e"'' $j_i(-s/i)) , 

(using Vj = 0,j < 0). Taking the Z-transform of this expression gives Y{(,,sh) = hB{(,,sh) V{£,)/{1 — e''^^) 
when ^ e^*'', where 

B{^,sh) = M-sh)+Y,[^A-sh) - e'''<f,-i{-sh)\ e • (3.9) 
It hence follows from (13.81) that 



and comparison with p.6p yields the alternative representation for the Z-transform of the weights qj : 

QiO = ^ fm(^)ds. (3.10) 



The expression B{£^, sh) /{I — e'''* ^) plays a role similar to that of {6{^)/h — s) ^ in standard CQ analysis, 
and it is the key quantity in determining whether the scheme is stable or not. Unfortunately it has a more 
complicated structure: it has an infinite vertical line of simple poles at s = Sk for /c G Z, where 

Sk ■■= I {- In ICI - * Arg(0 + i 27rfc) (3.11) 
h 

and the principal argument Arg(^) e (— 7r,7r]. Note that if |^| < 1 then Rc(sa:) > 0. 

To evaluate Q{(,) defined by (|3.10p for a given kernel function K{t) we can use either the left or right "D" 
contours illustrated in Figure [21 taking the limit i? — > cx) and setting 7 = lim/f^oo 7_r- Using the right contour 
gives 

QiO-±J<is.)BiC,s,h)-lun^±-^ j^T<is) (^) ds. 

The integral round does not necessarily vanish as i? — > cx) since for some basis functions (including higher 

order B-splines) the quantity — , = Oie'^"^) as Re(s) — > 00 for c > 1. We may also use the left contour 

1 — e"" 4 

when K{s) has simple poles at s = Kj with Rc(Kj) < and obtain the analogous result 

Q(0 = ^y^^^ lim(.-K,)ir(.)- hm A / Ki^,) ds. 

The asymptotic behaviour of the integral as i? — >■ 00 is determined primarily by K{s). The extension of 
this left contour approach to poles with higher multiplicity is straightforward. We illustrate the use of these 
formulae in Sec. 14.21 for various kernels K when the basis functions arc B-splines. 



A lm(s) 



X 




X 



Figure 2: The left and right "D" contours of radius R used for the stabihty and Z-transform calculations. The 
crosses are the poles p. lip and the vertical line of length (approximately) 2R is 77?. 

4 B-spline basis functions for (11.21) 

We now illustrate the theoretical framework introduced in Section |3] for basis functions which are B-splines 
on [0, oo). We begin by listing some general properties of B-splines which are needed in the subsequent analysis, 
and then examine the stability of the convolution spline approximation of ()1.2p for different example kernels. 
We also prove that the approximation given by p.3|) converges to the solution u of (|1.2p for general smooth 
a and K. The convergence rate is at most second order, no matter how high the polynomial degree, because 
quasi-interpolation by the Schoenbcrg B-splinc operator is at most 0{h^) (De Boor 1978). However a simple 
modification of the B-spline basis near t ~ can give higher order stable approximations of (II. 2[) . and this is 
analysed for the cubic case in Section [5l 



4.1 Notation and properties 

Throughout this section we assume that the basis functions {0j}j>o ^-rc (iso-geometric) B-splincs of polynomial 
degree m based on the nodes (or knots) tj for j > 0. We assume uniform spacing, so tj^i — tj ^ h for j > 0. 
It is necessary for the B-spline basis functions to have the sum to unity property ()2.12p in the whole interval 
[0, oo), and we introduce m new knots tj — for j — ~m : —1. The mth degree B-splincs are bj^{t) for 
j > —m, and B-splincs of degree m > are recursively defined in terms of those of lower degree as follows, 
using the convention that &™(t) = for j < — m. 

Definition 4.1. (De Boor 1978) When m = 

1 ifte [tj , tj+i) for j >0 and 



^^-^ 1 otherwise 



If m > then 



where the convention is that 0/0 is interpreted as 0. 

We shall make use of several B-splinc properties in Section 2] (see for example standard references such as De 
Boor (1978), Schumakcr (2007)), which we list here for convenience. 

B-spline properties 

PI. Compact support. b"^{t) = outwith [tj,tj+m+i), and h"^{tj) = unless j = — m. 



P2. Translate property. If j > then = ¥^{t/h — j), where the functions 6™ are defined recursively: 

771 + 1 — T 



b\T) = {l andifm>l: 6'"(t) ^ - 6'"-^(t) + ^ ^ ' _ i) . 

[0 otherwise ~ to ™ 

It follows that (/)j(t) = &™(r + 771 — j) for j > to. 

oo 

P3. Sum to unity. ^ = 1 for all t > 0. 



P4. Moments. 

ft. 



^ J + 1 ^j — m 
TO + 1 



and 



tb"l^it)dt 



(m + 1)(to + 2) 



m+l 



(jYl -|- T j (tTI ~\~ J ~]~ 1) 

P5. Shocnberg quasi-interpolation. Suppose that m > 1 and set t"^ — for 7 = —fn : —1 

■' 2m 

00 

and = ij+(m+i)/2 for j > 0. Then ^ t™ = t when t > 0. 

It follows from properties IP1[ IP3I and IP5I above that 

oo 

/W= E fitT)bT(*) + '^('^') (4-1) 

for any / e C^[0,oo), and if / £ C''+^[0,oo) for p > 2 and t e [te,ti+i) for some ^ > 0, then 



m- E /(^D^rw^E-^^'^^''^ 



fe=2 



Throughout this section we shall use basis functions 



and it follows from IP II that the CQ weights are 



3~m 



(t) forj>0, 



+ 0{hP+'^). (4.2) 



(4.3) 



(4.4) 



The convergence analysis relies crucially on knowing the values of the weights when K is a constant, and this 
follows immediately from IP4I when K =1 the weights qj of (|4.4p are given by 

^ ^ / t^T^ for j : TO - 1 
h \ 1 if J > TO . 

4.2 Stability results for convolution B-splines 

We now use the theoretical framework introduced in Section [3] to examine the stability of the convolution 
B-spline approximation of ()1.2p for different example kernels which capture some of the important properties 
of TDBIE problems. These are: Kit) equal to a constant, a step function, and the highly oscillatory kernels 
K{t) = Jo{u!t) or cos(aji), where to can be of the order of 1/h. We use B„i{£^, sh) to denote the function defined 
by (|3.9p for the degree to basis functions, and QmiO to denote the coefficient Z-transform given by p.lOp . 
The first few values of B„i{^, sh) are listed in Table [2] those for higher values of to are more complicated, but 
are easily computed in a standard algebraic manipulation package. 

In three of the cases QmiO is a rational polynomial in ^ and Theorem 13.11 can be used to determine stability. 
The Bessel function case is more complicated, and stability is determined from the Z-transform inversion 
formula by bounding the coefficients p„ of p.4p directly. Note that this bound is independent of n, and so is 
a practically useful stability result, in contrast with the (essentially) uncheckable hypotheses needed in Davies 
& Duncan (2003). 



V< Qn 1 1 n r> 

degree 


B (f 


1im r, R 'i) 


f> (p-s \ 
-Dm 1*7 


m = 


s-^(e" - 1) 


1 


3-^(6" - 1) 


m = 1 




(1+0/2 


s-'^ie" - l)2e-" 


TO = 2 


s--' {e^'ie -0 + e'(2 - 2C's + ^{s' + 2s - 2))+ 
((2s + 3)^-5^ -25-2-^2)) 




s-3(e« - l)3e-2^ 


m — 3 









Table 2: The function Bm{^,s) for m = : 3. See text for details. 



4.2.1 Constant kernel: K{t) = 1, transform /\(s) = 1/s 
Integrating p.lOp round the left contour in Figure [2] gives 

s^o 1 — e*"-^ (m + l)(l — 4)^ 

The function Qm has to simple roots on the unit circle, and stability of the approximation then follows from 
Theorem 13. II (Note that stability also follows from the convergence result of Sec. 14.31 ) 



4.2.2 Discontinuous step-function kernel: K{t) — 1 for t E [0,L], otherwise 0. 

Discontinuous kernels can arise in TDBIE problems, even when the scattering surface F is smooth and closed. 
Examples (in Laplace transformed representation) are given in Banjai & Lubich (2011, Sec. 6.1) and Sauter 
& Veit (20116, Sec. 4.1) describing time domain scattering where only the zeroth order harmonic in space is 
excited on the surface of a sphere. Similar, but more complicated discontinuous kernels are described in Sauter 
& Veit (2011 &) for more general scattering from spheres involving higher spatial harmonics. 
We assume that the duration L is independent of h and denote the integer part of L/h by M, i.e. when h is 
sufficiently small, L = {M + r) h for integer M > m and r S [0,1). It is simplest to work with the explicit 
Z-transform formula p.Sp using the weights given in (|4.4p . Results for to = : 3 are summarised below. 

Case m = 0: =re'+ ^" = 7^ i''^^'''' + - ^) - l) • 

n—Q 

When r E (0, 1) it can be shown that the M roots of Qq satisfy > 1 for j = 1 : M, and when r ~ 
there are M — 1 simple roots = exp(i27rj/A/) for j ^ 1 : M — 1. Hence Theorem 13.11 implies that the to = 
scheme is stable for all L. 



Case m = 1: 



r2 ^A^+i + (2 r - r2 + 1) ^'^^ + 2 J2n=i + 1 ior r e (0, 1) . 



h 



When r = the roots of Qi are —1 and exp{i2TTj /M) for j = 1 : A/ — 1, so there is a double root when M is 
even, giving a linear (and not exponential) instability. When r > the modulus of the product of the Af + 1 
roots of Qi{^) is l/r^ and it can be shown that there is a real negative root with modulus strictly greater than 
l/r^, implying that one or more of the roots has < 1 and the scheme is unstable. 
In summary, the to = 1 scheme is unstable unless L/h is an odd integer, which is not always possible. 

Case m = 2 : 3. Similar arguments to those above can be used to show that the scheme is unstable for any 
r in these cases. Although note that the modified cubic spline basis functions described in Section [5] give 
completely stable results for this kernel. 



4.2.3 K{t) = Jo{ujt), transform K{s) = l/^/s'^+cu'^ 

This is the kernel function that arises when considering TDBIE scattering from the flat surface M.'^, where w 
can be of the order of l//i (i.e. huj is bounded as /i — > 0, but does not necessarily tend to zero). Its Laplace 
transform has a branch cut between the values s — ±iui, and the Z-transform Qm(C) of weights is not a 



rational polynomial. We can still establish stability directly for the impulse response sequence {pn} defined in 
p.4p using a change of variable in the Z-transform inversion formula (Doetsch 1971, eq. 37.7) to get 



2tt 



,iny 



(4.6) 



where we have set ^ = e ^ with s ~ a + ii] and ct > and then changed to scaled variables x = ah and 
y ~ rjh. This yields the bound 



\qa\dy 



(4.7) 



\Q^ie—^y)\ 

when t„ < T, which holds for any fixed <t > when the singularities of the integrand arc to the left of x. 
Note that this bound is independent of ?i, and the scheme is stable at a given frequency a; if the integral term 
in (|4.6[) remains bounded as ft. — 0. This can be demonstrated using the right contour in Fig. [2] to calculate 
Qm(e^^'*) but it is more straightforward to work directly with p.7|) . 
It follows from standard properties of the B-splinc basis (j4.3|) that 



90 



where functions K'^ '^^ (<) are recursively defined by 
Note that 



K'^-^\t')dt' forfc^O, 1,.. 



for all t > 0. 

Properties of the B-spline basis functions can also be exploited to write (|3.7p as 

(1-0"+^ 



Z{K 



(-m-l) 



(4.8) 



(4.9) 



where the "correction" terms are Co = 0, Ci — 0, 



/l2 



2ft3 



The presence of these terms is because for m = : 1 the basis functions are pure translates, while for m > 2, 
there are different shaped basis functions at the start. The function K^^''^ has Laplace transform 

1 



if(-'=)(s) 



and it follows from the Poisson sum formula relating Z and Laplace transforms that 



where Sj = s + iiirj /h, and we use this expression in (|4.9p in order to bound the integral term in (|4.7p . 
When m = it is possible to obtain an analytic bound when uj < 7r//i, and a careful numerical approximation 
of the integral (|4.7p indicates that the p„ are bounded for w up to (at least) 20 7r/ft,. The situation is more 
complicated for m = 1 : 3, and in these cases we give numerical bounds. 



Case m = 0: From (l4?8l) and K9 



and 



ml 




IQo(e- 



1+ E ^(A"//o°; 

kei/o 



11 X]fe6z/fc 




1/2 



> l/o"l 



1+ E ^(/°//o; 

fcez/0 



It can be shown that 



nnn J] {fjQ = 1 + ^ ^ (/°//o) l.=o > ^ 
fcez/o fcGZ/o 



when < X < 1, < w/i < tt and |?;| < tt. In this case 



\qo\ ^ 3 Vx^ + 27r2 + y2 ^ 3 \/PT2^ Tre-^/^ 



|Qo(e-^-»a)| - 2 |e^+'''!' - 1| 



using Jordan's inequahty. Together with (|4.7p this proves that the scheme is stable in the sense of Definition 
l3.Il for frequency a; in the contiguous interval < ujh < tt. Numerical evaluation of the right hand side of (j4.7p 
indicates that the bound is 

\Pn\ < 1.3 e'^^ 

when h is sufBciently small (so that x < 0.1) and < uh < 207r. Further numerical tests computing p„ directly 
from (|3.4p for a finite number of steps n < 2500 and the same range of values of uh indicate that |p„| < 1, 
consistent with the estimate above. There is no indication of instability at any value of Loh tested and we 
speculate that this scheme is stable for all uj. 

Case m = 1 : Finding an explicit bound for the integral in (|4.7p is significantly more complicated and perhaps 
even intractible here so we only consider its direct numerical evaluation over a range of frequencies and values 
of X = ha close to 0. However there is an extra complication because qo/Qi{z) has a pole at z = —1. This 
is most obvious when we set cj = and get qo/Qi{z) = (1 — ^)/(l + z) from (|4.5p . When uj ^ there is no 
simple formula, but it is still possible to show by direct evaluation of the summation formula for (3i(e^^^*^) 
that the pole remains when < ujh < tt. The pole renders the bound in (|4.7p less useful since 



/I m^^\ ^ ^(^°S(1/^)) ^ oo as X ^ 



\Qie---'«)\ 

(where x = ha) and hence |p„| < (Co + Ci log(l/ft-))e°'^ as ft, ^ 0, which does not satisfy the stability 
requirement of Def. 13.11 Fortunately the singularity can be removed by writing 

10 a A , ,. 90(1 + 

— ' AP(4), where a = hm 



QiiO (1 + 1 QiiO 

so that AP(0 is bounded as ^ — )- —1. The sequence {p„} can then be written as p„ = a (—1)" + Ap„ where 
Apn is bounded in the same way as (|4.7p : 



|Ap„| <^J \AP{x + iy)\dy. 



(4.10) 



Numerical evaluation of the integral over frequencies < uth < tt indicates that |Ap„| < 1.1 and < a < 2 for 
nh < T and < x < 1/10. Combining this with direct evaluation of (|4.7p when ojh G [tt, 207r] and there is not 
a pole indicates that 

K|<Ce-^ where C = ( ^'J ' "J^S^'!), , 
'■^ ' - 1 1.1 , wft e [tt, 207r] 



for nh < T and < a; < 1/10 satisfying the stability Def. 13.11 Further numerical tests computing p„ directly 
from (j3.4p for a finite number of steps n < 2500 and the same range of values of ujh indicate that |pn| < 2 for 
ojh e [0,0. 77r) and |p„| < 1 for uih £ [0.77r, 207r], consistent with the estimate above. Again we speculate that 
this scheme is stable for all lu. 

Case m = 2: The function qo/Q2{0 appears to have two poles on the unit circle when ujh e [0, L) where L w 
2.55, symmetrically located at ^ = ^±tK'^h) ^ simple case u ^ 0, ((45|) gives Q2{£,) = qo (l+^+C^)/(l-C) , 

and so fj,(0) = 27r/3. Numerical evidence indicates that fi(ujh) > uih and that ^ increases until the two poles 
meet where I-l{L) = tt. At that point stability in the sense of Def. 13. II breaks down since there does not appear 
to be any compensating factor in the numerator to reduce the order of this double singularity. 
We locate the poles numerically, and remove them from the integrand qo/Q2{C) in a similar way to the previous 
case. The simplest form that captures the main features of the behaviour is 

10 «(1-C) 



uo e-'^icosfi+i 



maxJpJ,K(t) = jQ(cot) 



max^ Ip^l , K(t) = cos(co t) 



10' 



10" 




Figure 3: Plots of max{|p„| : < n < 2500} against ujh for the B-spline schemes with m = : 3 apphed to the 
highly oscillatory kernels K{t) = Jo{ujt) (left plot) and K{t) = coscot (right plot). See text for more details. 



so that by direct inversion of the Z transform 

Pn = a{cos{nfi) — sin(?i/.i) tan(/i/2)) + Ap„. 
For < w/k L R:! 2.55 we find that < a < 1 and from (|i?T(I)) that |Ap„| < 0.8, giving 

\Pn \ < 0.8 + sec(/i(a;/i)/2) 

for nh < T when < x < 1/10. This satisfies Def. 13.11 since 27r/3 < p,{u!h) < tt, but since sec(/^/2) — > oo 
as /X — >■ TT, the possibility for instability is clear. Further numerical tests computing pn directly from (|3.4p for 
a finite number of steps show very close and consistent agreement with this bound on \pn\, with instability 
appearing as predicted at ujh = L w 2.55 - i.e. there is a contiguous interval of stability ujh G [0, L) with 
L w 2.55. 

Case m = 3: Not surprisingly this case is more complicated still. When w = the three poles of Qo/QsiO are 
on the unit circle at ^ = —l,e^'''"^'^. However, when ujh > increases, the real-valued pole at ^ = —1 moves 
(harmlessly) outside the unit circle while the other complex conjugate pair moves inside causing instability. 
Numerical tests computing p„ directly from p.4[) for fixed vlaues of ujh show behaviour consistent with this: 
we see apparent stability for larger values of h which disappears as — >■ 0, i.e. this scheme is stable only when 
u is fixed (so that ujh — > 0). 

Similar results can be proved for the (more straightforward) oscillatory kernel K{t) = cosut, as summarised 
below. 

• m = : the scheme is stable at any frequency uj for which uh £ [0, tt) ; 

• m = 1 : the scheme is stable at any frequency uj for which cj/i G [0, 27r) ; 

• m = 2 : the scheme is stable at any frequency lu for which ujh G [0, 9), where 6 = 1.9747 . . . ; 

• m = 3 : there is no 0(1) interval of stability for ujh, but the scheme is stable for bounded oj. 

The stability results for highly oscillatory kernels are illustrated in Figure [3] The plots show max„ |p„ | for 
n = : 2500 for the B-spline schemes with m = : 3 applied to the kernels K{t) = Jo{ojt) (left plot) and 
K{t) = cosujt (right plot). Over the range ut G [0,7r/ft,] shown, the general stability behaviour for these two 
kernels is similar. In particular the left plot illustrates the stability when m = 0, 1, while scheme m = 2 is 
stable for uih G [0, L) with L « 2.55. On the right plot, scheme m = is stable except at ujh ~ tt, scheme 
m = 1 is stable and m ~ 2 scheme is stable for uth G [0, L) with L w 1.97. On both plots the 771 = 3 scheme is 
clearly unstable when lj = 0{l/h). 



4.3 Convergence results for ( 13. Sh 



Formal convergence of a method when applied to a smooth VIE problem is certainly a necessary condition for 
it to behave well (for VIEs or TDBIEs) , and we now prove that the collocation spline approximation p.3p with 
B-splinc basis functions (|4.3p converges to the solution u of (|1.2p for general smooth a and K. The analysis 
proceeds by considering the case K = 1 and then using Taylor expansion to show that the same result also 
holds for smooth K with K{0) = 1 when h is small enough (see e.g. Brunner (2004)). Thus it does not apply 
to the important case of an oscillatory kernel where the oscillation frequency uj — 0{l/h), whose stability was 
analysed above for the Bcssel function and cosine kernels. 

When m = the approximation p.2p is the same as using piecewisc constant collocation (at the interval 
endpoints), and this has been fully analysed (see e.g. Brunner (2004) for details). Here we assume that 
m > 1 (note that this includes the well-known case of piecewise linear approximations of (|1.2p ). and show that 
convergence is always second order, no matter how high the polynomial degree, because quasi-interpolation by 
the Schoenberg B-splinc operator is at most 0{h?) (De Boor 1978). This is in marked contrast to discontinuous 
polynomial collocation or Galcrkin approximations of (|1.2p which converge at optimal order (Brunner 2004, 
Brunner, Davies & Duncan 2009, Brunner, Davics & Duncan 2012). However a simple modification of the 
B-spline basis near t = Q can give higher order stable approximations of (|1.2p . as illustrated when to = 3 in 
Sec. El 

The approximation error e„(t) for n > 0, i > is 



en{t) = U{tn -t)-Y. ' (4-11) 

J=0 



where the coefficients Vj satisfy 



a{tn) ^^qjVn~j (4.12) 
i=o 

for weights qj as defined in (|4.4p . Note that = (because a(0) = 0) and so the sums above can be taken 
from J = to n — 1, and it then follows from Property [FT] that e„(t) = for t > t„ and each weight can be 
written as qj = /q" K{t) bj^_^{t) dt . Hence, multiplying (14.111) by K{t) and integrating gives 

/ K{t)en{t)dt^ K{t)u{tn-t)dt-y^qjVn-3=Q: 

Jo Jo ^-^0 

by (|1.2p and (|4.12p . i.e. e„ is orthogonal to K on {0,t„). The formal convergence result is as follows. 
Theorem 4.1. Suppose that m > 1 and the conditions lll.5\) and h2.1\) hold for d > 4. Then 

\en{t)\<Ch^ (4.13) 
for t G [tyn^T], for some C independent of n and h. If m = 1 then {^.IS^ holds for t G [0,T]. 
Remarks: 

• The m — 1 case has been fully analysed (Brunner 2004) and is just included for completeness. 

• The restriction to second order convergence for m > 1 is a fundamental aspect of quasi-interpolation by 
classical B-splines and not an artefact of the proof, and is illustrated in Figure H] when (|2.ip holds with 
d = 4. 

• Equation (|4.13p trivially holds for t > t„ (because e„(t) = 0), and so it is enough to prove the result for 
t € [tm,tn) when TO > 1, where n<T/h. 



Proof. We first express e„(t) in terms of coefficients ;= u(tfc_|_(,„_i)/2) ~ ^fc • Substituting the quasi- 
interpolation result (|4.ip with f{t) = u(<„ — t) in the definition (|4.1ip of e„ {t) gives 

n n+rn 
j=m j=n+l 




Figure 4: Convergence results for the approximation of ()1.2p for two smooth kernel functions with maximum 
time r = 10 and a{t) = t^exp(— 50(t — 1/2)^). Convergence rates of 0{h^) for splines of degree m > 1 and 
0{h) for m = are clear. Stability results for the highly oscillatory kernels coswt and Jo(wi) where the 
frequency lu can be 0{l/h) are given in Section 14.2.31 and Figure [3l 



where we have used &J1,„(0 = for j < m and j > n + m. It follows from the assumptions (|1.5p and (|2.ip 
that u{ch) = 0{h'^) for any constant c. This implies that the second sum term in the previous equation is 
0{h'^), and hence yields 

n 
j=m 

for t G [tm, tn) with tn < T . Because there are at most m + 1 nonzero terms in this sum for any i, it is sufficient 
to show that there exists a constant C independent of h such that 

|ej| < Ch^ for aU j < T/h. (4.14) 



To prove (|4.14p note that it follows from (14.121) that qj Vn^j ~ / u(i„ — t) dt and so 

Yl T ^"-J ^ II f "(*n-J + (™-l)/2) - T / ^(0 ^itn ~ t) dt Rn . (4.15) 

If h is sufficiently small, then expanding K{t) and using lP4l gives 

1^-1 + hK'{0) + 0{h') for j=0: m - 1 

~r ~ ) m + 1 2(to + 1)(to + 2) (4.ibj 

[ K{tj^rn) + ^h{m+l)K'{t,^rn)+0{h^) if J > m . 

It then follows from the quasi-interpolation result ()4.2p with p ~ 3 that when n > m, 

ri-l 3 , -.Nfcii »t.,, I £+m 



hRn^Y^T.^' (<„-,)+0(/i*) for 77,^ = / X(t) <^ {t - Uf - (^^-™ - ^^)' ™(^) \ ■ 

£=0 fe=2 ' "^'f I j=£ ' 



It is then straightforward to show Rn+i — 2 i?„ + i?„_i = C'(/i'^) and after some manipulation using (|4.16p the 
second central difference of (|4.15p can be written as 

rn n 
e=0 i=m+l 

where 7„ ~ 0{h'^) and all the bltc bounded. This can be written as a one-step recurrence for the vector 
gn g Rm+i -vvith components (5" = £„+j for j — : m. The recurrence is 5"'^^ — for j = : m — 1 with 



S'^-'^^ given by (|4.17p . which gives the matrix- vector system 



n—l — m 

6!'+^ ^ {M + hWo) S!" + Wn^iT^^ + lne"' (4.18) 

where e™ = (0, . . . , 0, 1)^, each matrix is bounded and M £ R(™+i)^(™+i) is the circulant matrix whose 
only nonzero components are Mm.o = Mjj+i — 1 for j — Q : m — 1. The eigenvalues of M are the (to + 1)— th 
roots of unity, and are hence distinct. Following Brunner (1978) we note that M belongs to Ortega's (1990, 
§1.3) Class M, and so there is a vector norm |j • ||, on M™+^ for which the induced matrix norm satisfies 
= p{M) = 1. Taking this norm of (|4.18p then implies that there is a constant C such that 

Ti+l — m 

and the top bound of (|4.14p gives < C h'^. Standard arguments can then be used to show that ||^"||* < f„ 

where = C h'^ and 

n 

Vn+i^{l + Ch)vn + Ch^Yvi + Ch^ for7z>0. (4.19) 

1=0 

This has the solution = A+X\ + A_ Xi where A± = 1 + 0(h) and A± = Oih"^). Hence < Ci e'^^ ^ 
for n < T/h, which concludes the proof of (|4.14p and hence (|4.13p . □ 



5 Modified cubic spline basis functions for (11.21) 

The results of the previous section illustrate that the convolution spline framework can be used to derive new 
VIE approximations and prove their stability in cases not covered by standard convergence analysis (such 
as discontinuous or highly oscillatory kernels), but the restriction to second order convergence for the B- 
spline basis (|4.3p is not competitive with RK-based CQ methods (Banjai 2010, Banjai & Lubich 2011, Banjai 
et al. 2011). We now show how a slight modification of the to = 3 B-spline basis near t = can yield methods 
which are more accurate and have better stability properties than (|4.3p . 

5.1 Derivation 

It is simpler to define the modified basis functions in terms of B-splines centred at zero, and we set B{t) — 
b^{t + 2), so supp(i3) = (—2,2). For j > 3 the basis functions are (j)j{t) = B{t — j), and we choose (pj for 
j —Q -.2 to satisfy the "parabolic runout" conditions at t = 0, namely 

0o(i) =S(t)+3B(t + l), (t>i{t) = B{t-l)-'iB{t + l), (j)2{t) = B{t - 2) + B{t + I) . 

This means that 

oo 

^ /((/)j(t)- B(t -.?•)) = for < e [0,oo) for r = : 2, (5.1) 

where wc set 0_i(i) = 0, and in particular it ensures that quasi-interpolation in terms of {4'j}JLo linearity- 
preserving on [0,oo). Weights qj for the VIE (|1.2p are given in terms of these basis functions by (|2.10p and 
the approximation at t = i„ is 

n 
3=0 

where the Vj coefficients are given by p.3p . The individual coefficients can be directly obtained from Un{t — tn) 
by introducing dual basis functions $fc(0 such that 




There are many ways in which a suitable dual basis can be chosen, and one possibility is to use continuous 
piecewise cubic functions on [0,oo) defined with respect to knots in N. For fc > 3, we set ^k{t) = ^{t — k) 



where $ is the even continuous piecewise cubic function on [—2, 2] which is zero at ±2 and is at the interior 
knots ±1, and satisfies 

$(i - £) B{t) dt = Soj for ^ = : 3, 

'-2 

and it is also possible to find suitable continuous piecewise cubics $j, with support in [0, fc + 2) for A: = : 2. 
Calculating these dual basis functions in an algebraic manipulation package is straightforward, although their 
coefficients are messy and we omit their details. Once the $fc have been obtained it can be shown by Taylor 
expanding that for any sufficiently smooth function /, 

k+2 ^2 

fitn - hs) ds = f{tn-k) - -T /"(^n-fc) + 0{h^) . (5.3) 

max(0,fc-2) O 

Multiplying (|5.2p by ^^{t/h) and integrating over [0,oo) gives 

rk+2 

Vn-k = / Unitn - hs) ^k{s) ds 

Jmax(0,fc-2) 

and we now use this representation and (|5.3p in order to obtain the best approximation of the solution u of 
(|1.2p in terms of the basis functions 0j when (|1.5p holds with d = A. Specifically, we set 

n 

<tn -t)=Y^ Un-j (j)j{t/h) + Rn{t) , (5.4) 

where 

w(i„ - hs) $j(s) ds = u{tn-j) - — u"(i„-j) + 0{h'^) . (5.5) 



max(0,i-2) 6 



We now show that this scheme is fourth order accurate, and then discuss its stability properties for discontin- 
uous and highly oscillatory kernels. 

5.2 Convergence 

It follows from (|5.4[) - (j5.5|) that the approximation error e„(t) := u{tn — t) — Un{tn — i) for this scheme is 

n 

en{t) = Y,^n-j 4>j{tlh) + R^it) 

3=0 

where Sk = Uk — Vk- We now show that modifying the basis functions near t = as described above improves 
the scheme's accuracy from second to fourth order. 

Theorem 5.1. Suppose that the conditions U.5]) and h2.1\l hold for d > 6. Then there exists a constant C 
independent of n and h such that if tn < T, the VIE approximation error satisfies 

\enit)\<Ch^ (5.6) 

for t e [ti,i„]. 

Proof. It is straightforward to verify that if s G [0, 1), 

2 

Rnitk + hs) = u{tn-k - hs) - ^ Un^k~l <i>k+l{k + s) 

where 0-1 is taken to be zero, and hence it follows from (|5.ip . (|5.5p and standard B-spline properties that 

~ r h\\-sf^'\t^)l%^0{h^) whenfc^O 

i?„(tfc + M-| 0(;^4) ^henfc>l. 

It is thus sufficient to show that each = 0(h^) (because at most four basis functions are nonzero for any t). 
The proof follows that of Thm. 14.11 and the expression analogous to (|4.15p is 



n— 1 n „i 

j=0 fe=o-^o 



Taking the second central difference, using the fact that when K = 1 the scheme coefficients are 



go = 5/1/8, (?i= 5/1/6, 92= 25/1/24, and qj = h for j > 3 (5.7) 

gives an expression like (|4.18p with m = 3 where the bottom row of the matrix AI is now [—1, 6, 0, 10]/15. The 
eigenvalues of M are all distinct and its spectral radius is p{M) = 1, and this again yields a bound Vn which 
satisfies a difference scheme like (j4.19|) whose final term is C . This gives an 0{h'^) bound for each \ek\ with 
tk < T, and ([51]) follows. □ 



5.3 Stability for discontinuous and highly oscillatory kernels 

We first consider the discontinuous kernel introduced in Section [4.2.21 The duration L — [M + r)h is fixed 
independent of h such that integer M > 5 and re [0,1). In this case the qj arc given by (|5.7p for j = : M — 2, 
the coefficients qj for j = M — 1 : M + 2 are polynomial in r and satisfy 

qM-2 = h > qM-i > qM > qM+i > qM+2 > , 

and qj = for j > M + 3. 

The stability proof follows that of (Davies & Duncan 2013, §3.1.3); forward differencing p.4[) gives 

A/+3 , _ X 

15p„ + 5p„_i + 5p„_2-p„-3 + 24 J2 = forn>2, 

where the first stability factors are po = 1, Pi = ^Qi/lo = —4/3 and we set pj = for j < 0. Similar 
manipulations to those in the previous subsection then give 

11 8 

\Pn\< Zn-1 + g Z„_(A/_i) for 71 > 2 

where now = maxj<k \Pj \ ■ It then follows from similar arguments to those used in (Davies & Duncan 2013, 
§3.1.3) that if n > M - 1 then 

i.e. Zn < 6 Zn-{M-i)y and |p„| < Z„ ~ 4/3 for 71 = 1 : A/ — 2. In combination these give the stability bound 

for all n < T/h when h is sufficiently small, and so the stability coefficients are bounded independently of h 
as required by Def. 13.11 Note that this bound is very pessimistic, and in practice the factor is about 1.5^/^. 
This scheme is also stable when applied to the oscillatory kernels Jo(wt) and cos(a;t) examined in Section l4.2.3l 
The analysis is much simpler to carry out for this scheme, since there is no problem with poles of 1/Q{C) on 
or near the unit circle |^| = 1. We find that 



\Pn\ < Ce"^ where C = | J'^ ' 



2.2 , K{t) = cos(wi) 
K{t) = Joiiut) 



for all ujh S [0, 207r] and beyond. This is verified in tests computing p„ directly from (|3.4p for a finite number of 
steps n < 2500 and the same range of values oiuih. They indicate that \pn\ < 1.82 (for cos(aji)) and \pn\ < 4/3 
(for J(){ujt)), consistent with the estimate above. 



5.4 Numerical results 

Numerical comparisons of the new modified cubic B-spline approximation for (|1.2p are compared with the 
CQ BDF2 method and with the convolution B-splincs with m = 0, 3 from Sec. S] in Figure [S] when a satisfes 
(|2.ip with d = A. The convergence rates for a smooth kernel (left subplot) are as expected: CQ BDF2 has 
rate 0{h^); B-splines and 3 have rates 0{h) and 0{h?)\ and the modified B-splinc 3 has rate 0{h'^). The 
discontinuous kernel problem (right subplot) has the same convergence rates, apart from the isogeometric 
m — 3 approximation which was shown to be unstable in Section [4. 2. 2| and is not illustrated. 




Figure 5; Convergence results for the approximation of (jl.2[) with maximum time T = 10 and a(t) = 
t^exp(— 50(t — 1/2)^) (see text for details). The stability and 0{h'^) convergence rate of modified B-spline 3 
is clear for both problems. 

6 Convolution splines for TDBIEs 

The results shown here are obtained by approximating the solution of the TDBIE in space by piccewise 
constant basis functions on a generally irregular triangular grid, and in time by convolution splines as 

n AI 

U{x,tn - t') ~^^ulr^ 4>3{t' Ih) Vk{x) , 

where the arc the modified cubic B-splincs from Sec. 15. 1[ 

, , J 1 for a; e Tfe 
Vk[x)-<^ otherwise 

and Ffc is the fcth triangle on the surface F. The spatial Galcrkin formulation of the problem gives the 
time-marching scheme (|1.3|) with 



Qrn ff ^rai\x y\/h) ^ ^^^^ j ^^^^^n^^^^ 

JJr.xr^ \x-y\ Jr, 

The matrices Q"* are symmetric and calculating each component involves a four dimensional integral. The ofi- 
diagonal components and the components a" have smooth integrands and are approximated by a composite 
triangular quadrature with 16 sub-trianglcs, each of which is sixth order with 12 quadrature points. Each 
diagonal element of the Q™ has a singular integrand, and is first converted into smooth subintegrals using a 
Duffy-type transformation, and then approximated by the same quadrature rule as the rest of the calculation. 
Similar results are obtained when piecewise linear spatial basis functions are used. 

Figure [6] shows results for (jl.ip when F is a flat plate and a sphere, both with incident field a{x,t) ~ ao{t + 
to — |a;|)/|a;|, where ao{t) — exp(— 20(t — 1/2)^) for t > 0. In these two tests the mesh ratio is chosen to be 
h/Ax = 1/2 , where Ax is the size of a typical space mesh element. The left-hand graph shows the size (and 
hence stability) of the approximate solution on a unit square plate. The growth in the Loo norm is due to a 
corner singularity (Holm, Maischak & Stephan 1996) in the exact solution while the maximum of the 1-norm 
is well-behaved as the mesh is refined. The right-hand plot shows the maximum relative error (i.e. the error 
normalised by the maximum size of the solution) for scattering from a unit sphere - an exact solution for this 
problem is given in Sauter & Veit (20116). The dotted lines show a second order convergence slope, and this 
is the best which can be expected from the spatial approximation, despite the higher order accuracy of the 
temporal approximation. 

Figure [7] demonstrates the impact of changing the mesh ratio h/ Ax in the sphere scattering example described 
above. The left plot shows the maximum relative error in the solution against the number of space elements, 
and the right plot shows the dependence of this error on h. As one would expect for a scheme with higher 
order accuracy in time than space, the time error decreases faster than the space error when the mesh is 
refined with fixed mesh ratio, and the space error eventually dominates. This is clear on the left plot for 
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Figure 6: Solution plots for the TBIE p.ip when the scattering surface F is a unit square plate (left plot) and 
a unit sphere (right plot). See text for details of the space and time approximations used. 




Figure 7: Solution plots for the TBIE for mesh ratios /i/Aa; ranging from 0.25 to 4 when the scattering 
surface F is a unit sphere. The dotted line indicates second order convergence. See text for details of the space 
and time approximations used. 



ratios h/Ax = 0.25,0.5 and 1. For the larger mesh ratios the time step size is simply too big to resolve the 
input function accurately, and the asymptotic convergence regime has not yet been reached. Increasing the 
mesh ratio decreases the number of time steps (and hence matrices Q'") used, but increases the number of 
non-zero entries in each matrix. However the net result is a decrease in the computational cost for both the 
time marching calculation and the computation of the matrices Q'". 



7 Conclusions 

We have derived a new framework for time-stepping approximations of the TDBIE (|l.ip . The system matrices 
Q™ of the resulting scheme (|1.3p have the same degree of sparsity as for a space-time Galerkin approximation, 
but are much more straightforward to calculate, especially when higher order (smoother) B-splines are used. 
The method is constructed as an approximation scheme for the VIE p.2p , and key properties are its backward 
time aspect (j2.9p . that the basis functions have compact support and are (mainly) translates, and they satisfy 
the sum to unity condition (|2.12p . 

These properties permit a full stability analysis for VIEs with kernels which capture some of the important 
properties of TDBIE problems, as illustrated for B-spline basis functions in Secs. l4.2 land l5.3l In particular the 
analysis gives a stability bound for the pn coefficients for the Besscl function kernel (|1.4|) when a; w 1/h. This 
means that the convergence proof in Davies & Duncan (2003) for the TDBIE (jl.ip on F = could be applied 
to an approximation which is a Fourier interpolant in space and a convolution spline in time, without the need 
to impose an additional (essentially uncheckable) stability assumption. The modified cubic convolution spline 



approximation of Sec. [S] is fourth order accurate and gives a very stable approximation for discontinuous and 
highly oscillatory VIE kernels. The TDBIE numerical test results indicate that the scheme is very stable and 
performs well. 

There is current interest in TDBIE time-stepping methods which can use variable time-steps. See e.g. (Lopez- 
Fernandez & Sauter 2011, Lopez-Fernandez & Sauter 2012) for convolution quadrature methods and (Glafke 
2013) for space and time adaptation in the full space-time Galerkin method for scattering problems in 2D 
space. Because our TDBIE time-stepping method is based on B-splines (whose key approximation properties 
are retained for a non-uniform knot distribution) and standard piecewise polynomials in space, the various 
strategies described in (Glafke 2013) for space and time adaptation could be applied. However, using variable 
time-stepping in any TDBIE approximation algorithm imposes an overhead because it essentially involves 
recalculating the Q™ matrices of (|1.3p at every time step, and this is extremely expensive. 
We note that there are also other choices of basis functions which seem to give stable approximations of (|1.2p 
and (|l.ip when used in the same convolution framework. Non-polynomial temporal basis functions are 
introduced in Davies & Duncan (2013); they are translates for j > 2, and for j = : 1 they arc defined as 
described in Wu & Schaback (1994) for radial basis function (RBF) multi-quadrics in order to ensure that quasi- 
interpolation in terms of {4'j}JLo linearity-preserving. They also work well as a temporal approximation of 
the TDBIE (|l.ip . but because the basis functions are global the system matrices need to be sparsified (but 
this is straightforward because they are highly peaked). The method derived in (Davies & Duncan 2013) is 
second order accurate, and the fourth order modified cubic B-spline approximation of Sec. [5] is a significant 
improvement. Extension of this approach to modified B-splincs with m > 3 is work in progress. 
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